subroutine set_ic( dfu, uxMat, uyMat, f0vec, finit, it)    ! outputs

  use dimensions
  use mpi_wrap

  use fft_wrap, only   : fft_forward, fft_backward
  use fft_3d

  use parameters, only : kz, vPhase, ux0, uy0, &                      ! real,     phys
                         nz, iWaveFrame, iModx, iMody, iModz, &       ! integer,  phys
                         ModAmp,  &                                   ! real,     phys
                         restart, &                                   ! integer,  num
                         pi, iunit                                    ! constants

  implicit none

!-- interface input/output --

  real,    intent(out) :: dfu(jv, jz, jy, jx, jDir)
  real,    intent(out) :: uxMat(jz, jy, jx, jDir)
  real,    intent(out) :: uyMat(jz, jy, jx, jDir)
  real,    intent(out) :: f0vec(jv)
  complex, intent(out) :: finit(jv,jz,jDir)
  integer, intent(out) :: it


!-- local variables ---

  integer              :: ix, iy, iz, iv, iDir
  real                 :: MoDinZ(jz), MoDinX(jx), MoDinY(jy) 
  complex              :: dummy3D(jz,jy,jx)


!-- create vector "f0vec" --


  f0vec(:) = ( 1/sqrt(2*pi) )*exp(-0.5*(vvec(:) + vPhase*iWaveFrame)**2)
                 


!-- read restart or create new --

  if (restart) then

       write(msg,*) char(10), ' Reading restart', char(10)
       call mpi_wrap_msg
       call mpi_wrap_msg

       call read_bin(dfu, uxMat, uyMat, it)        


  else
 
       write(msg,*) char(10), ' Creating initial data', char(10)
       call mpi_wrap_msg

       it = 0
       dfu=0.0

       do iDir = 1, jDir
          uxMat(:,:,:,iDir) = ux0(iDir)
          uyMat(:,:,:,iDir) = uy0(iDir)
       enddo
 
       if ( iModx ) then
          MoDinX(:) = cos( dkx * xvec(:) )
       else
          MoDinX(:) = 0
       end if

       if ( iMody ) then
          MoDinY(:) = cos( dky * yvec(:) )
       else
          MoDinY(:) = 0
       end if

       if ( iModz ) then         
          MoDinZ(:) = cos(  kz * zvec(:) )
       else
          MoDinZ(:) = 0
       end if


       do ix=1,jx                                                      
          do iy=1,jy
             do iz=1,jz
                do iv=1,jv
                   dfu(iv,iz,iy,ix,:) =  &
                        ModAmp * f0vec(iv) * ( MoDinZ(iz) + MoDinX(ix) + MoDinY(iy))
                enddo
             enddo
          enddo
       enddo

  end if

!-- save initial (1D, equilibrium) dfu.

  do iDir = 1,jDir
     do iv=1,jv

        dummy3D(:,:,:) = dfu(iv,:,:,:,iDir)

        call fft_forward_x(dummy3D)
        call fft_forward_y(dummy3D)

        fInit(iv,:,iDir) = dummy3D(:,1,1)

     end do
  end do



  !fInit(:,:,:)=dfu(:, :, 1, 1, :)      

end subroutine set_ic
